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Abstract. Using a high performance computer cluster, we run simulations regarding 
an open problem about d-dimensional critical branching random walks in a random IID 
environment The environment is given by the rule that at every site independently, with 
probability p > 0, there is a cookie, completely suppressing the branching of any particle 
located there. 

The simulations suggest self averaging: the asymptotic survival probability in n 
steps is the same in the annealed and the quenched case; it is where q := 1 — p. 
This particular asymptotics indicates a non-trivial phenomenon: the tail of the survival 
probability (both in the annealed and the quenched case) is the same as in the case of 
non-spatial unit time critical branching, where the branching rule is modified: branching 
only takes place with probability q for every particle at every iteration. 
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1. Introduction 

In [3] a spatial branching model has been studied, where the underlying motion is a 
d-dimensional (c? > 1) Brownian motion, the particles perform dyadic branching, and the 
branching rate is affected by a random collection of reproduction suppressing sets dubbed 
mild obstacles. In fact the obstacle configuration was given by the union of balls with fixed 
radius, where the centers of the balls form a Poisson point process. The radius r plays no role 
in the results, but the Poisson intensity v > does. The main result of |3] is the quenched 
Law of Large Numbers for the population for all d > 1. The environment w (with law P'^) 
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has the property that the branching rate is f3i inside the obstacles and ^2 outside of them, 
with < Pi < ^2- The branching process given the environment uj was denoted by Z.(uj) 
and the total population size at i > was denoted by \Zt{uj)\. Define the average growth 
rate by rt — rt{u}) := l2S.^ifc£li, [3J it was shown that for almost every environment, 

lim {\ogtf/'\rt - P2) = -c{d, v) P'^ -probability, 

where c{d, v)>0 is an explicit constant. (This is a kind of LLN, because the expectation of 
the total population size obeys the same asymptotics.) 

It has also been shown that the branching Brownian motion with mild obstacles spreads 
less quickly than ordinary branching Brownian motion, and an upper estimate for its radial 
speed has been provided. 

More general offspring distributions (beyond the dyadic one considered in the main the- 
orems) were also discussed in liSj. In particular, the following question was posed. Consider 
the model where the offspring distribution is critical. One can easily prove (see Theorem 
12.21 below) that, despite of the obstacles, the system still dies out with probability one. 

Problem 1.1. What is the rate of decay for the survival probability? Is it still of order 
C/n as in the obstacle-free (non-spatial) case? 

In the present paper we are going to investigate this problem in a discretized setting. 
More precisely, we consider a modified version of the model, by replacing the Poisson point 
process with IID probabilities on Z**, as follows. Given an environment, the initial particle, 
located at the origin, first moves according to a nearest neighbor simple random walk, and 
immediately afterwards, the following happens to her: 

(1) If there is no cookie at the new location, the particle either vanishes or splits into 
two offspring particles, with equal probabilities. 

(2) If there is a cookie at the new location, nothing happens to the particle. 

The new generation then follows the same rule in the next unit time interval and produces 
the third generation, etc. 

Let p G [0, 1]. In the sequel Pp will denote the law of the cookies and P'^ will denote 
the law of the BRW given the environment uj. So, if Pp denotes the 'mixed' law in the 
environment with cookie probability p, we have 

Pp(.)=EpP-(-). 

Following standard terminology. P'^ and Pp will be called quenched and annealed probabil- 
ities, respectively. 

We close this section with a remark which is essentially taken from ^ with slight adap- 
tations. 

Remark 1.2. An alternative view on our setting is as follows. Let K denote the (random) set 
of those lattice points where there is a cookie. Then, our model can be viewed as a catalytic 
BRW as well — the catalytic set is then K"^ (in the sense that branching is 'made possible' 
there). Catalytic spatial branching (mostly for superprocesses though) has been the subject 
of vigorous research in the last twenty years initiated by Dawson, Fleischmann and others 
— see the survey papers [6J and [2J and references therein. In those models the individual 
branching rates of particles moving in space depend on the amount of contact between the 
particle ('reactant') and a certain random medium called the catalyst. The random medium 
is usually assumed to be a 'thin' random set (that could even be just one point) or another 
superprocess. Sometimes 'mutually' or even 'cyclically' catalytic branching is considered 

m- 

Our model is simpler than most catalytic models as our catalytic/blocking areas are fixed, 
whereas in several catalytic models they are moving. On the other hand, while for catalytic 
settings studied so far results were mostly only qualitative we are aiming to get quite sharp 
quantitative results. 
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For the discrete setting there are much less results available. One exampl43 is [S] , where 
the branching particle system on TJ^ is so that its branching is catalyzed by another au- 
tonomous particle system on TL'^ . There are two types of particles, the A-particles ('catalyst') 
and the B-particles ('reactant'). They move, branch and interact in the following way. Let 
Af^(a;, s) and Nb{x, s) denote the number of A- (resp. B-)particles at a; G Z'^ and at time 
s G [0, oo). (All Na{x, 0) and Nb{x, 0), x G Z"^ are independent Poisson variables with mean 
fiA (mb)-) Every yl-particle (i3-particle) performs independently a continuous-time random 
walk with jump rate Da (Db)- In addition a i3-particle dies at rate 5, and, when present 
at X at time s, it splits into two in the next ds time with probability I3Na{x, s)ds + o{ds). 
Conditionally on the system of the ^-particles, the jumps, deaths and splitting of the B- 
particles are independent. For large /3 the existence of a critical 6 is shown separating local 
extinction regime from local survival regime. 

Finally, note that while the continuous equivalent of an IID trap configuration on the 
lattice is a Poisson trap configuration on M.'^, there is an important difference between the 
two. The discrete setting has the advantage that the difference between the sets K and K'^ 
is no longer relevant. Indeed, in the discrete case the complement is also IID with a different 
parameter (self-duality), whereas in the continuous setting this nice duality is lost as the 
'Swiss cheese' K'^ is not the same type of geometric object as K; the latter is the case, for 
example, in [3]. 

2. Preliminaries 

In this section we present two simple statements which are relatively easy to verify rig- 
orously. Let Sn denote the event of survival for n > 0. That is, Sn ~ {Zn > 1}, where Z„ 
is the population size at time n. 

Theorem 2.1 (Monotonicity). Let < p < p < 1 and fix n > 0. Then 

Pp('S'n) 5: Pp(>S'n). 

Proof. First notice that it suffices to prove the following statement: 

Assume that we are given an environment with some 'red' cookies and some 
additional 'blue ' cookies. Then the probability of Sn with the additional 
cookies is larger than or equal to the probability without them. 

Indeed, one can argue by coupling as follows. Let q := 1 — p, 6 :^ p ^ p. First let us 
consider the cookies that are coming with IID probabilities and parameter p. These will 
be the 'red' cookies. Now with probability 6/q at each site independently, add a blue 
cookie. Then the probability for any given site, that there is at least one cookie there is 
p + 5/q~pS/q~p + S ~ p. Now delete those blue cookies where there was a red cookie too. 
This way, the red cookies plus the additional blue cookies together correspond to parameter 
P- 

We are thus going to prove the statement in italics now, using an argument due to 
S. Kuznetsov. Consider the generating functions of no branching and critical branching: 
(pi{z) — z and <f2{z) = ^(1 -I- z'^), and note that ipi < ip2 on R. Fix an environment and 
define 

u{n,x,N) := Pn^xiS^), 

that is, the probability that the population emanating from a single particle which is at time 
n is located in x, becomes extinct at time N. Then, if the particle moves to the random 
location $„+i, one has 



A further example of the discrete setting is [T]. 
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where pi{^n+i) is the probabihtjQ of producing i offspring at the location ^n+i, and ip^"+^ 
is either ipi or ip2. 

Now consider two environments: one with only red cookies and another one where there 
are also some additional blue cookies, and let us denote the corresponding functions by ui 
and U2- We have 

ui{n, X, N) ^ Eip\-+' [uiin + 1, £,n+i,N)] 

and 

U2{n,x,N) = E^l-+' [u2 (n + 1 , C„+i , N)] . 
Clearly ui{N, x, N) — U2{N, x, N) = 0. Therefore, using that ipi < ip2 and by (backward) 
induction, U2 > ui for all n = 0, 1, — 1. In particular, 

Ui{0,X,N) < M2(0,X,iV), 

finishing the proof. □ 

We now give a precise statement and a rigorous proof concerning the result about eventual 
extinction mentioned briefly above. 

Theorem 2.2 (Extinction). Let < p < 1 and let A denote the event that the population 
survives forever. Then, for ¥p- almost every environment, P'^{A) — 0. 

Proof. Let again Z„ denote the total population size at time n for n > 1. Then Z is a 
martingale with respect to the canonical filtration {J-"„;n > 1}. To see this, note that just 
like in the p — case, E{Zn+i — Zn \ J-'n) = 0, as the particles that do not branch (due to the 
presence of cookies) do not contribute to the increment. Being a nonnegative martingale, 
Z converges a.s. to a limit Zoo, and of course Z^c is nonnegative integer valued. We now 
show that for Pp-almost every environment, P^{Zoo = 0) = 1. Introduce the events 

• Ak {Zoo - fc} for fc > 1, 

• B: branching occurs at infinitely many times < cri < (T2 < •■• 

Clearly, A = Ufe>iAfe = {Zoo > !}• We first show that 

(2.1) for Pp - a.e. environment, P'^{B''A) = 0. 

Clearly, it is enough to show that 'Pp{B'^A) — 0. 

Now, B'^A C C, where C denotes the event that there exists a first time N such that 
for n > N, there is no branching and particles survive and stay in the region of cookies. 
On C, one can pick randomly a particle starting at N, and follow her path; this path visits 
infinitely many points P'^-a.s. whatever w ifl Since this path is independent of w and p < I, 
the Pp-probability that it contains a cookie at each of its sites is zero. Hence Pp(C) = 0, 
and (P?T|) follows. 

On the other hand, for each fc > 1, there is a < 1, such that the probability that 
the population size remains unchanged (it remains fc) at am is not more than pk for every 
m > 1, uniformly in lo. Thus, 

P"(BAfe) = P"(B n = fc for all large enough to}) = 0, 

whatever uj is. Using this along with (|2.ip . we have that for Pp-almost every to, P'^{Ak) = 
P'^iB'^Ak) + P'^iBAk) =0, fc > 1, and so P"(A) =0. □ 



So either pi = 1 or po = P2 = 1/2, according to whether there is no cookie there or there is one. 
■^Because for every oj, it is true P"-a.s., that every particle that does not branch, has a path that visits 
infinitely many points 
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Algorithm 1 The annealed simulation function of the code. 



typedef struct { 

Tcell cell; 
int iter ; 
} Tparticle ; 

typedef vector < Tparticle > Tparticles 



location 

number of iterations 



// dimension 
// cookie prob ability 

// maximum number of iteration steps 
// longest survival 



void simulate ( 
int dim , 
double p , 

int maxiter , 
int &iter 
) { 

board . clear {) ; 
iter = ; 

Tparticles particles ; 

particles . reserve ( maxiter / 10) ; 

Tparticle particle ; 

particle . cell=Tcell {dim , 0); 

particle . iter = 0; 

particles . push_ back {particle ) ; 

while { \ particles . empty () iter < maxiter) { 

int coord = mtr.randint {dim — 1); 
int dir = 2 * mtr.randint (1) — 1; 
particles . back (). cell [coord] 
if (cookie (p, particles . back 
p art i cl es . back {).iter++; 

iter = max { iter , particles . back {) . iter 

} 

else 

if {mtr2 . randint (1)) { 

particles . back {).iter++; 
iter = max (iter, particles . back ().ite 
particles . push_back { p art i cl e s . back () 

} 

else 

particles . pop back (); 

} 

} 



+= dir ; 
O.ceJJ)) { 



remove cookies 
iteration counter 
all the particles 
/ reserve room 

cell at origin 

particle at center 
alive particles 
random coordinate 
-1 or +1 
move 

there is a cookie 
survived iteration 
longest survival 

no cookie 
prob ability 0.5 
survived iteration 
longest survival 
duplicate particle 

die 

remove the particle 



3. Implementation 

The code for our simulations are written in Ch — h using the MPIqueue parallel library 
[5]. We ran the code on 96 cores using a computing cluster containing Quad-Core AMD 
Opteron(tm) 2350 CPU's. We used an implementation [T^ of the Mersenne Twister [7J to 
generate random numbers. The total running time for the simulations was several months. 

3.1. Annealed simulation. Algorithm 1 shows the C++ function that runs a single an- 
nealed simulation. We are essentially implementing a 'depth- first search.' Below we give a 
detailed description of the code. 

• line 1: We define a data type to store particles. 

• line 2: The location of the particle is stored in the cell field, that is a vector with 
the appropriate dimensions. 

• line 3: The iter field stores the number of iterations survived by the particle. 

• line 5: We define a data type to store all the living particles. 

• line 7: The simulation function takes three input variables and one output variable. 

• line 8: The dimension of the space is the first input. 

• line 9: The probability of a cookie at any given location is the second input. 

• line 10: The maximum number of allowed iterations is the third input. 
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• line 1 1 : The output of the function is the maximum number of iterations any particle 
survived. 

• line 13: We erase all the cookies from the board. Every run of the simulation uses 
a new cookie placement. 

• line 14: The initial value of the output must be zero. 

• line 15: We define a variable to store all our alive particles. 

• line 16: We reserve some space to store the particles. Making the reserved space too 
small results in unnecessary reallocation of the variable which degrades performance. 
On the other hand, reserving too much space can be a problem too since different 
CPU's compete with each other for RAM. 

• lines 18-19: The initial particle starts at the origin before the iterations start. 

• line 20: At the beginning we only have the initial particle. 

• line 21: We run the simulation while we have alive particles and none of them stayed 
alive for the maximum allowed number of iterations. 

• lines 22-23: We generate a random direction. 

• line 24: We move the last of our alive particles in this random direction. 

• line 25: We call the cookie function to check if there is a cookie at the new location of 
the particle. The cookie function checks in the global variable board if any particle 
already visited this location and as a result wc know already whether there is a 
cookie there. If no particle visited this location before, then the function uses the 
cookie probability to decide whether to place a cookie there or leave the location 
empty. This information is then stored for future visitors. 

• line 26: If there is a cookie at the new location, then the particle survived one more 
iteration, so we increment the iter variable. 

• line 27: It is possible that this is the longest surviving particle so far, so we update 
the output variable. 

• line 29: If there is no cookie at the new location, then the particle splits or dies. 

• line 30: Wc generate a random number to decide what happens. 

• lines 31-32: If the particle splits, then it survives so we update information about 
the number of iterations. 

• line 33: The particle splits, so we place a copy of it into our collection of particles 
as the last particle. 

• lines 35-36: If the particle dies, then we remove it from our collection of particles. 

The rest of our code takes care of the parallelization, data collection and the calculation 
of survival probabilities. The program splits the available nodes into a boss node and 
several worker nodes. The boss assigns simulation jobs to the workers. The workers call 
the simulation function several times. The boss node collects the results of these jobs and 
calculates the survival probabilities using all the available simulation runs. More precisely, 
Pp('S'n) is estimated as the number of simulation runs with longest survival value not smaller 
than n divided by the total number of runs. 

3.2. Quenched simulation. The code for quenched simulation is essentially the same with 
only minor modifications. In this version, line 13 of he simulation function is missing, since 

we do not want to replace the board at every simulation. 

The other change in the simulation function is at line 25. In the annealed case, every 
worker node has a local version of the board and the cookie function can create the board 
on the fly. In the quenched case, the worker nodes need to use the same board, so the 
cookie function cannot generate the board locally. The new version of the cookie function 
still stores information about the already visited locations. On the other hand, if a location 
is not visited yet, then the worker node asks the boss node whether this new location has 
a cookie. The boss node first checks whether the location was visited by any other particle 
at any other worker node. If the location was visited, then the boss already has a record 
of this location. Otherwise, the boss node uses the cookie probability to decide whether 
the location should have a cookie. Essentially, the boss node has the ultimate information 
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Figure 4.1. Results for an annealed 2-dimensional simulation. The graph 
shows the reciprocal of the survival probability as a function of the number 
of iterations. Each line represents a different cookie probability. One such 
line is the result of 10* runs of the simulation with a newly generated cookie 
landscape. 

about the board, but the worker nodes keep partial versions of the board and only consult 
the boss node when it is necessary. 

Remark 3.1. In the quenched case, note that if p„ denotes the relative frequency of survivals 
(up to n) after r runs for a fixed environment lu, that is. 



then using our method of simulation, the random variables p„ and pm are not independent 
for n ^ m, because the data are coming from the same r runs. 

Similarly, in the annealed case, for a fixed environment and a fixed run, the random 
variables Island 15^^ are not independent for n ^ m. 



We ran our annealed simulation on Z'' with d G {1, 2, 3}. The 1-dimensional case turned 
out to be the most challenging. So we start the description of our results with the 2- 
dimensional case. The 3-dimensional case produced essentially the same output as the 
2-dimensional case. 

4.1. Annealed simulation on 7?. We executed 10* runs allowing a maximum of n,nax = 
10000 iterations with p G {0, 0.025, 0.05, 0.075, 0.1, 0.2, . . . , 0.9, 0.925, 0.95, 0.975}. For p = 1 
we used the known survival probability of 1. Preliminary results made it clear that the 
simulation is more sensitive for small and large values of p. This is why we picked more of 
these values instead of a uniformly placed set of values. The reciprocal of the calculated 
survival probabilities are shown in Figure |4T1 The figure suggests that n n> l/Pj,(S'„) is 
asymptotically linear. We calculated the slopes for these curves from the values at 4ninax/5 
and nniax- These slope values are shown in Figure l4?2] 
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Figure 4.2. Results for an annealed 2-diniensional simulation. The graph 
shows the apparent slopes in Figure |4?11 (i.e. the limits of the functions of 
Figure 15. ip as a function of the cookie probability together with the graph 
oip^ (1 -p)/2. 
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Table 1. Exact and simulated survival probability values Pp(S'„) on I?. 
The simulated values are calculated from the data used in Figure 14.11 



To verify the correctness of our simulation we computed the exact theoretical survival 
probabilities after the first two iterations. It is easy to see that Pp(5'i) = 1/2 + p/2 and 

Pp(^2) = 3/8 + llp/32 + 3p^/m + 3pV32. 
Table [T] compares some of the exact and simulated values. 

4.2. Annealed simulation on Z^. A 1-dimensional simulation with 10^ runs and n^ax = 
10000 produces less satisfactory results as shown in Figure l473l The reasons behind this are 
explained in Subsection 15.21 below, where a discussion is given concerning the fluctuations 
of the empirical curves in the figures. Essentially, in the annealed case, small values of 
Pp{Sn) result in large errors (see Subsection 15.21 for more explanation) and therefore we 
modified the original algorithm by introducing a stopping rule: when the estimated value of 
Pp('S'n) reaches a certain small threshold value, we stop and do not simulate more iterations. 
Fortunately, when larger threshold values needed, they are actually large: we obtained slower 
convergence for large values of p, and, clearly, for those values, the probability Pp(5„) is 
large. The threshold value was set 1/4000, based on trial and error. This way, we stopped 
the iteration at nstop(p); the slopes were then calculated from the values at 4nstop(p)/5 and 
nstop(p)- See Figure IIjI 

After adjusting the algorithm by using the above stopping rule, the curve indeed straight- 
ened out and the picture became very similar to the 2-dimensional one in Figure l4?2l 
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Figure 4.3. Results for an annealed 1-dimensional simulation. Every pa- 
rameter for this simulation is chosen to be the same as that of Figure 14.21 
except the dimension. 
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Figure 4.4. Annealed 1-dimensional simulation with 959965800 runs. For 
small values of p (lines at the top) , a small iteration number would actu- 
ally give better results, because otherwise the survival probability Pp(5„) 
becomes too small, even with a huge number of runs. On the other hand, 
for large p values (lines at the bottom) one needs large iteration numbers 
because the convergence is apparently slow. The squares represent the it- 
eration thresholds after which p„ < -mKn- 
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Table 2. Exact and simulated survival probability values Pp(S'„) on Z 
The simulated values are calculated from the data used in Figure 14.41 
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Figure 4.5. Results for a quenched 2-dimensional simulation with 10^ 
runs. Each line represents a different cookie landscape. One such line 
is the result of 10^ runs of the simulation. The lines are in three groups 
corresponding to three different cookie probability. Each group has 50 lines. 
The cookie probabilities from top to bottom are 0.25, 0.5 and 0.75. The 
total number of simulations required for this graph is 3 • 50 • 10^. The total 
running time was about 29 hours. 



To verify the correctness of our simulation we computed the exact theoretical survival 
probabilities after the first two iterations. It is easy to see that Pp{Si) = 1/2 + p/2 and 

Pp{S2) = 3/8 + 5p/16 + pV4 + p^lG. 

Table [3] compares some of the exact and simulated values. 

4.3. Quenched simulation. From the annealed simulation it has been clear that conver- 
gence is much faster in two dimensions than in one dimension. Therefore, in the quenched 
case we chose to present our results in the d = 2 case. In fact, we got qualitatively similar 
results for d = 1. 

In Figure 14.51 we see three bundles corresponding to three values of p. Those bundles are 
very thin, so essentially the same thing happens for every realization, and the slopes of the 
lines are roughly 3/8, 1/4 and 1/8 from top to bottom, corresponding to p = 0.25, p — 0.5 
and p = 0.75, respectively. That is, for each one of these values of p, the slope is the same 
as in the annealed case. 
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Although Figure [1751 is about the d = 2 case, we have a similar simulation result for d = 1; 
in fact we conjecture that this qualitative behavior (that is, the coincidence of the first order 
asymptotics of the quenched and annealed survival probability) will hold for all d > 1 . 

5. Interpretation of the simulation results 

5.1. Main finding. Recall first the classic result due to Kolmogorov [4, Formula 10.8], that 
for critical unit time branching with generating function 95, as n — 00, 

2 

Pfsurvival up to n) ^ 

n(f"[l) 

As a particular case, let us consider now a non-spatial toy model as follows. Suppose that 
branching occurs with probability q G (0, 1), and then it is critical binary, that is, consider 
the generating function 

ip{z) = {l~q)z + ^q{l + z^). 
It then follows that, as n — > 00, 

2 

(5.1) P(survival up to n) ^ — . 

qn 

Turning back to our spatial model, the simulations suggest (Figures 14.11 and l 4.5p the self 
averaging property of the model: as explained in the previous section, the asymptotics for 
the annealed and the quenched case are the same. In fact, this asymptotics is the same as 
the one in 115.1]) . where p = 1 — q is the probability that a site has a cookie. In other words, 
despite our model being spatial, in an asymptotic sense, the parameter q simply plays the 
role of the branching probability of the above non-spatial toy model. To put it yet another 
way, q only introduces a 'time-change.' 

The intuitive picture behind this asymptotics is that there is nothing either the environ- 
ment or the BRW could do to increase the chance of survival, at least as far as the leading 
order term is concerned (as opposed to well known models, for example when a single Brow- 
nian motion is placed into random medium [9J). Hence, given any environment the particles 
move freely and experience branching at q proportion of the time elapsed (quenched case), 
and the asymptotics agrees with the one obtained in the non-spatial setting as in (|5.ip . 
Furthermore, creating a 'good environment' (annealed case) and staying in the part of the 
lattice with cookies for very long would be 'too expensive.' 

Note that whenever the total population size reduces to one, the probability of that 
particle staying in the region of cookies is known to be much less than 0{l/n) (hard obstacle 
problem for random walk) . So the optimal strategy for this particle to survive is obviously 
not to try to stay completely in that region and thus avoid branching. Rather, survival will 
mostly be possible because of the potentially large family tree stemming from that particle. 
In fact, the formula Pp(S'„) ^ together with the martingale property of \Zn\, implies 
that Ep(|Z„| I Sn) -^-n. 

Notice that the straight lines on Figures [4.21 and [4.31 start at the value 1/2, that is, as 
p 4- 0, one gets the well known non-spatial asymptotics 2/n, which is a particular case of 
Kolmogorov's result above. We conclude that there is apparently no discontinuity at p = 
(no cookies) for the quantity lim„_).oo nP(survival up to n). 

5.2. Interpretation of the fluctuations in the diagrams. Since we estimated the recip- 
rocal of the survival probabilities and not the probabilities themselves both in the annealed 
and the quenched case (Figures 14.11 and [ 4 . 5 p . one cannot expect good approximation results 
when those probabilities are small. Indeed, in the annealed case, if /9„:=Pp(S'„) (with p 
being fixed) and pn denotes the relative frequency obtained from simulations, then the Law 
of Large Numbers (LLN) only says, that if the number of runs is large, then the difference 
\pn —pn\ is small. However, looking at the difference of the reciprocals 
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Figure 5.1. Results for an annealed 2-dimensional simulation with 10® 
runs. The graph shows the reciprocal of the survival probability divided 
by the number of iterations as a function of the number of iterations. The 
data used to create the graph is the same as that of Figure 14.11 



it is clear that a small p„ value magnifies the error; in fact the effect is squared as p„ is close 
to pm exactly because of the LLN. This effect is the reason of the 'zigzagging' of the line 
on Figure for small p values. In fact, small p values result in small /?„ values in light of 
Theorem 12.11 and the continuity property mentioned at the end of the previous subsection. 
Clearly, there is a competition between p„ being small (as a result of p being small and n 
being large) on the one hand, and the large number of runs on the other. The first makes 
the denominator small in , while the second makes the numerator small according to 

LLN. 

Looking at Figure one notices another peculiarity in the 1-dimensional setting. For 
large values of p, the empirical curve is slightly under the straight line. The explanation 
for the relatively poor fit is simply that the iteration number is not large enough for the 
asymptotics to 'kick in.' 

These arguments are bolstered by the experimental findings that increasing the number 
of runs helps for small p values, whereas increasing the number of iterations helps for large 
ones. For example, in Figure[43]we increased the maximal iteration number rimax to 200000 
and plotted n H> (nPp(S'„))^^. One can see that for small p values it is beneficial to stop 
the iterations earlier, but for large values large iteration numbers give better results. 

We do not have an explanation, however, for the deviation downward from the straight 
line (for large p values) in Figure We suggest, as an open problem, to find at least a 
heuristic explanation for this phenomenon. 

Interestingly, for higher dimensions there is apparently a perfect fit for large values of 
p, indicating that for higher dimensions the asymptotics is much quicker then for d — 1. 
Figure [5T] checks the assumption (for d — 2, annealed) that the reciprocal of the survival 
probability is ^ + o(n) as n — oo, by dividing the reciprocal of the survival probability by 
n. The graphs convincingly show the existence of a limit, depending on p. 
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Figure 5.2. Zooming in at the bottom part of Figure [5TT] (i.e. large p 
values): the convergence is apparently from below. 
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Figure 6.1. Annealed 1-dimensional simulation with 7,259,965,800 runs 
and p = 0.5. 
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Figure 6.2. Annealed simulation with p = 0.5. The solid curve shows the 
l-dimensional result, the dashed curve shows the 2-dimensional result. 



6. Beyond the first order asymptotics 

In this section we will attempt to draw conclusions about more delicate phenomena 
beyond the first order asymptotics, and those conclusions will necessarily be less reliable 
than the ones in the previous sections. 

6.1. Two dimensions. Consider again Figure [SH] Zooming in gives Figure [5T2l Looking 
at Figure [44l Figure [Ql and Figure [521 foi' small p values (top lines) the convergence seems 
to be from above, and for large p values it seems to be from below. 

6.2. One dimension. For d = 1 we obtained figures somewhat similar to the 2-dimensional 
ones, which we summarize below without actually providing them. Simulation seems to 
suggest that for not too small values of p. the convergence is also from below; this is in line 
with the fact that, as we have already discussed, on Figure the l-dimensional empirical 
curve is below the straight line for large p values. For very small p's, the direction of the 
convergence is not clear from the pictures. Although the convergence is apparently quicker, 
the effects are 'blurred' due to the magnification of error explained earlier. 

The following conjecture is based on Figure [Ql 

Conjecture 6.1 (Second order asymptotics). The l-dimensional annealed survival proba- 
bility obeys the following second order asymptotics: 

Pp{Sn) = — + fin), 
nq 

where lim„_j.oo /('^) • '^^^^ = C > 0, and C may depend on p. 

6.3. Comparison between one and two dimensions. The annealed convergence to the 
limit 2/q seems to be quite different for d = 1 and d — 2. Figure [6?2l shows this difference, 
and in particular, it illustrates that in 1-dimension, the convergence is slower and it is 
apparently from below for p = 0.5. 



CRITICAL BRANCHING RANDOM WALK 



15 



References 

1. Sergio Albeverio and Leonid V. Bogachev, Branching random walk in a catalytic medium. I, Basic 
equations, Positivity 4 (2000), no. 1, 41-100. 

2. D. A. Dawson and K. Fleischmann, Catalytic and mutually catalytic super- Brownian motions, Seminar 
on Stochastic Analysis, Random Fields and Applications, III (Ascona, 1999), Progr. Probab., vol. 52, 
Birkhauser, Basel, 2002, pp. 89-110. 

3. Janos Englander, Quenched law of large numbers for branching Brownian motion in a random medium, 
Ann. Inst. Henri Poincare Probab. Stat. 44 (2008), no. 3, 490-518. 

4. Theodore E. Harris, The theory of branching processes, Dover Phoenix Editions, Dover Publications 
Inc., Mineola, NY, 2002, Corrected reprint of the 1963 original [Springer, Berlin; (29 #664)]. 

5. Harry Kesten and Vladas Sidoravicius, Branching random walk with catalysts. Electron. J. Probab. 8 
(2003), no. 5, 51 pp. (electronic). 

6. Achim Klenke, A review on spatial catalytic branching. Stochastic models (Ottawa, ON, 1998), CMS 
Conf. Proc, vol. 26, Amer. Math. Soc, Providence, RI, 2000, pp. 245-263. 

7. Makoto Matsumoto and Takuji Nishimura, Mersenne twister: a 623-dimensionally equidistributed uni- 
form pseudo-random number generator, ACM Trans. Model. Comput. Simul., no. 1, 3-30. 

8. John M. Neuberger, Nandor Sieben, and James W. Swift, MPIqueue: A simple library implementing 
parallel job queues in C-I-+, (preprint). 

9. Alain-Sol Sznitman, Brownian motion, obstacles and random media. Springer Monographs in Mathe- 
matics, Springer- Verlag, Berlin, 1998. 

10. Richard Wagner, http://www-personal.umich.edu/'^wagnerr/MersenneTwister.html. 

Department of Mathematics, University of Colorado, Boulder, CO-80309-0395, and De- 
partment OF Mathematics and Statistics, Northern Arizona University, Flagstaff, AZ-86011. 
E-mail address: Janos.Englander@Colorado.edu, nandor.sieben@nau.edu 

URL: http: / /euclid. Colorado . edu/~englaiidj/MyBoulderPage .html , |http : // j an . ucc . nau . edu/~ns46] 



